setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20210321.H3K分析")
group1=c("117012","13649","8544","807","200")
group2=c("Fetal_Brain_Female","Fetal_Brain_Male")
group3=c("enhancer","promoter","all","")
i=1
sel1=list.files(pattern = "Fetal_Brain",full.names = T,path = group1[i])
j=1
tmp=read.csv(sel1[j],header = T)

rt=tmp[-1,]
for (i in 1:length(group1)) {
  for (j in 1:length(group2)) {
    for (k in 1:length(group3)) {
      fn=paste0(group1[i],"/","consider.all.VAF.",group2[j],"DNase.all.peaks.",group3[k],"_statis.csv")
      tmp=read.csv(fn,header = T)
      tmp$Term=group1[i]
      tmp$Tissue=group2[j]
      tmp$regions=group3[k]
      rt=rbind(rt,tmp)
    }
  }
}

library(ggplot2)
library(ggsci)
library(multiplot)
i=1

p=list()
for(i in 1:length(group1)){
  for(j in 1:length(group2)){
    tmp=rt[rt$Term==group1[i],]
    tmp=tmp[tmp$Tissue==group2[1],]
    tmp$opposite.ratio=tmp$opposite/(tmp$same+tmp$opposite)
    tmp$group=ifelse(tmp$pvalue<0.05,"S","N.S")
    tmp1=data.frame(ratio=tmp$same_ratio,type="same.ratio",regions=tmp$regions,group=tmp$group)
    tmp2=data.frame(ratio=tmp$opposite.ratio,type="opposite.ratio",regions=tmp$regions,group=tmp$group)
    data=rbind(tmp1,tmp2)
    data$group=factor(data$group,levels = c("N.S","S"))
    data$regions=factor(data$regions,levels = c("enhancer","promoter","","all"))
    
    p[[2*(i-1)+j]]=ggplot(data=data,aes(x=type,y=ratio,group=regions,color=regions))+
      geom_line(size=1)+geom_point(aes(shape=group),size=3)+
      theme_light()+labs(x="",y="Ratio",title = paste0(group1[i]," ASH"))+
      theme(axis.title.y=element_text(size=16),axis.text.y=element_text(face="bold",size=15),axis.text.x=element_text(size=15))+
      scale_color_npg() 
  }
}

layout <- matrix(c(1,2,3,4,5,6,7,8,11,9,10,11), nrow = 3,byrow = T)
multiplot(plotlist=p,layout=layout)


ggplot(data=data,aes(x=type,y=ratio,group=regions,color=regions))+
  geom_line(size=1)+geom_point(aes(shape=group),size=3)+
  theme_light()+labs(x="",y="Ratio",title = paste0(group1[i]," ASH"))+
  theme(axis.title.y=element_text(size=16),axis.text.y=element_text(face="bold",size=15),axis.text.x=element_text(size=15))+
  scale_color_npg() 

rt.plot <- rt
rt.plot = rt.plot[rt.plot$pvalue<0.05,]
rt.plot[rt.plot$regions=="",]$regions = "Not enhancer or promoter"
rt.plot$opposite.ratio <- 1-rt.plot$same_ratio
rt.plot$Group <- paste0(rt.plot$Term, "_",rt.plot$Tissue)
rt.plot=rt.plot[!rt.plot$Term=="13649",]	##去掉甲基化的结果
rt.plot=data.frame(Same=rt.plot$same_ratio,Opposite=rt.plot$opposite.ratio,Group=rt.plot$Group,regions=rt.plot$regions,Tissue=rt.plot$Tissue)
rt.plot1=rt.plot[,-1]
colnames(rt.plot1)[1] <- "ratio"
rt.plot1$Ratio <- "Opposite"

rt.plot2 <- rt.plot[,-2]
colnames(rt.plot2)[1] <- "ratio"
rt.plot2$Ratio <- "Same"

str1 <- paste(group1,"Fetal_Brain_Female",sep="_")
str2<- paste(group1,"Fetal_Brain_Male",sep="_")

rt.plot <- rbind(rt.plot1, rt.plot2)
rt.plot$regions <- factor(rt.plot$regions, levels = c("enhancer", "promoter", "Not enhancer or promoter", "all"))
rt.plot$Ratio <- factor(rt.plot$Ratio, levels = c("Opposite", "Same"))
rt.plot$Group <- factor(rt.plot$Group, levels = c(str1,str2))

ggplot(data = rt.plot, 
       aes(x = Ratio, y = ratio, group = regions, color = regions))+
  geom_line(size=1)+geom_point(size=3)+
  #scale_shape_manual(values = c(1,2)) +
  theme_light()+labs(x="",y="Ratio",title = "")+
  theme(axis.title.y=element_text(size=16),axis.text.y=element_text(face="bold",size=15),axis.text.x=element_text(size=15))+
  scale_color_npg() +
  facet_wrap(~Group, scales = "free")

p <- paste0("p",1:20)
for (i in 1:length(rt.plot.Group)) {
  assign(p[i], ggplot(data = rt.plot[which(rt.plot$Group == rt.plot.Group[i]),], 
                      aes(x = Ratio, y = ratio, group = regions, color = regions))+
           geom_line(size=1)+geom_point(aes(shape=Sign),size=3)+
           scale_shape_manual(values = rt.plot[which(rt.plot$Group == rt.plot.Group[i]),"Sign1"]) +
           theme_light()+labs(x="",y="Ratio",title = paste0(rt.plot.Group[i]," ASH"))+
           theme(axis.title.y=element_text(size=16),axis.text.y=element_text(face="bold",size=15),axis.text.x=element_text(size=15))+
           scale_color_npg() )
}

cowplot::plot_grid(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10, ncol = 4)
